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Extreme mass ratio inspirals (EMRIs) (capture and inspiral of a compact stellar mass object 
into a Massive Black Hole (MBH)) are among the most interesting objects for the gravitational 
wave astronomy. It is a very challenging task to detect those sources with the accurate estimation 
parameters of binaries primarily due to a large number of the secondary maxima on the likelihood 
surface. Search algorithms based on the matched filtering require computation of the gravitational 
waveform hundreds of thousands of times, which is currently not feasible with the most accurate 
(faithful) models of EMRIs. Here we propose to use a phenomenological template family which 
covers a large range of EMRIs parameter space. We use these phenomenological templates to 
detect the signal in the simulated data and then, assuming a particular EMRI model, estimate 
the physical parameters of the binary. We have separated the detection problem, which is done 
in a model-independent way, from the parameter estimation. For the latter one, we need to adopt 
the model for inspiral in order to map phenomenological parameters onto the physical parameter 
characterizing EMRIs. 



I. INTRODUCTION 

Stellar compact objects like a black hole, neutron star 
or white dwarf in the cusp surrounding the massive black 
hole (MBH) in the galactic nuclei could be deployed on 
a very eccentric orbit due to N-body interaction. Such 
an object could either plunge (directly or after few or- 
bits) into MBH or form an EMRI: inspiralling compact 
object on originally very eccentric orbit which shrinks 
and circularizes due to loss of the energy and angular 
orbital momentum through gravitational radiation. The 
compact object spends ~ 10 6 orbits in the very strong 
field of a MBH before it plunges, all this orbital evolu- 
tion will be encoded in the phase of emitted gravitational 
waves (GWs). Space based GW observatories, like LISA 
or similar planned missions, will observe those sources 
few years before the plunge. By fitting precisely the GW 
phase one can extract extremely accurate parameters of 
a binary system [T] (like mass and spin of MBH M, a, 
mass of a small object to, inclination of the orbital plane 
(to the spin of MBH) , orbital eccentricity and semi-latus 
rectum (io> e O)Po) & t some fiducial moment of time to, 
location of the source on the sky (0, cf) and more). 

Precise tracking of the GW phase implies that we can 
also test the nature of the central massive object. The 
general belief is that it should be a MBH with surround- 
ing spacetime described by a Kerr solution. The nature 
of the spacetime affects the orbital evolution of the com- 
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pact object which in turn could be extracted from the 
GW phase. Kerr spacetime is described by only two 
parameters: black hole's mass and spin, as stated by a 
"no-hair" theorem. The spacetime could be decomposed 
in the multipole moments of a central massive object, 
and, for Kerr BH, all moments depend only on M and 
a: Mi + iSi = (ia) l M l+1 where M; and Si are mass and 
current moments. We could measure three first moments 
(mass, spin and quadrupole moment) [2 , and check the 
"Kerrness" of a spacetime. In general, the deviations 
from Kerr could come in several ways: (i) it is Kerr BH 
but there is an additional perturber (gas disk, another 
MBH) (ii) it is not Kerr BH but some other object sat- 
isfying GR (boson star, gravastar), (hi) there are devia- 
tions from GR. For discussion on the topics we refer the 
reader to (3H5] and references therein. 

Modeling orbital evolution even within GR is not yet 
fully complete. Large mass ratio allows us to consider a 
small compact object as a perturbation on the Kerr back- 
ground spacetime, and treat the problem perturbatively 
in orders of the mass ratio. In zero order approximation 
the compact object moves on a geodesic orbit, however, 
as soon as we assign the mass to it, it creates its own 
gravitational field interacting with the background and 
this system emits gravitational radiation. The force re- 
sulting from the interaction of the self field with the back- 
ground is called self force, and the motion of the compact 
object could be seen as the forced geodesic motion. Al- 
ternative interpretation is that the motion is governed 
by a geodesic motion but in the perturbed spacetime. 
Calculation of the self force is a complicated task which 
is accomplished for the orbits around Schwarzschild BH 
only [51 [7] , the Kerr spacetime is underway. There are 
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also questions concerning the calculation of the orbital 
evolution under the self force: the self force depends on 
the past history of the compact object (which is usually 
assumed to be a geodesic in the background spacetime). 
To compute the motion under the self force one can use 
the osculating elements approach [8J, or self-consistent 
approach of direct integration of the regularized equa- 
tions [5]. For more details on this subject we refer to 

ma. 

All in all, the modeling of the orbital evolution and 
the GW signal is a complex task which requires signif- 
icant theoretical and computational developments. The 
latter prevents us currently from using the state-of-art 
GW models of EMRIs in our data analysis explorations. 
In majority of the cases the phenomenological model sug- 
gested in p], so called "analytic kludge" (AK), is used. 
It is based on Post-Newtonian expressions and puts to- 
gether all relevant physics of EMRIs. However, this 
model has restrictions in the number of harmonics and in 
their strength, and any search algorithm which relies on 
its specific harmonic content will not work for a more re- 
alistic model of GW signal. The main motivation of this 
work is to create the phenomenological search template 
family which would fit a very large range of EMRI-like 
signals. The typical EMRI signal consists of a set of 
harmonics of three (slowly evolving) orbital frequencies, 
and we will use it as a basis of our template. The phe- 
nomenological template consists of harmonics with 
constant amplitude and slowly evolving phase which we 
decompose in a Taylor series. Truncation of the Taylor 
series and the assumption about constant amplitude set 
restrictions on the duration over which the phenomeno- 
logical template can fit an EMRI signal. The amplitude 
of EMRI's harmonics changes due to shrinking of the or- 
bit (overall amplitude increases), circularization of the 
orbit (power is shifted to lower harmonics) and slight 
change in the inclination of the orbit to the spin of MBH. 
Using more terms in the Taylor series helps to track phase 
of the EMRI signal for longer time (which is more im- 
portant than accurate description of the amplitude). Fi- 
nally, we decide on the number of harmonics to use in 
the template (and their indices) based on the analysis of 
the harmonic structure of the Numerical Kludge (NK) 
model llj of EMRI in different parts of the parameter 
space. The restriction that the phenomenological wave- 
form (PW) is valid only for a limited period of time is 
very weak since we can fit the signal piecewise, as long 
as the accumulated signal-to-noise ratio (SNR) over that 
time is significant to claim presence of the signal. In this 
work we consider only those parts of the EMRI signal 
where the orbital frequencies are not decreasing which is 
true over almost all time of the inspiral and breaks quite 
close to the plunge. However, this is not really necessary 
since we did not restrict the values of frequency deriva- 
tives to positive values during the search. 

The PW family is quite generic and does not depend on 
the orbital evolution, or, in other words, the orbital evo- 
lution of the binary is encoded in the Taylor coefficients 



of phase of each harmonic. This allows us to detect an 
EMRI signal in a model independent way. Once the har- 
monics of the signal are recovered we can analyze them 
using a specific EMRI model to recover physical parame- 
ters of the system. It is at this point we need the orbital 
evolution with high accuracy, which involves computa- 
tion of the self-force and tests of possible deviations from 
the "Kerrness". 

After constructing the phenomenological waveform we 
perform blind searches on the simulated data without 
noise (to avoid stochastic errors in the parameter estima- 
tion) and with the noise. We have used the NK waveform 
(as described in [TT]) as a model of our signal and the 
orbital evolution according to [12]. We have also used 
Markov chain Monte-Carlo (MCMC) search with phe- 
nomenological waveforms on the simulated three month 
of data. This search has provided us with multiple local 
maxima in the likelihood which we gathered and ana- 
lyzed in a similar way as described in |13j . We associate 
local maxima in the likelihood with partial detections of 
the signal and construct the time-frequency map of the 
detected (patchy) harmonics of the source. The next step 
is to assume the model for the orbital evolution and, by 
matching the found time- frequency tracks to the harmon- 
ics of the signal, estimate parameters of the binary sys- 
tem. We have used the same model for the orbital evolu- 
tion as in the simulated data sets and recovered physical 
parameters with precision better than few percent. 

The paper is organized as follows. In the next Section 
we will give a brief overview of available models for GWs 
from EMRIs. In Section |ITT[ we introduce PW family in 
details. We describe MCMC search with PWs in Sec- 
tion |IV| Analysis of MCMC results and mapping to the 
physical parameters are done in the Section [V] Finally 
we conclude with a summary Section [VI| 

II. REVIEW OF EMRI WAVEFORMS 

As was already mentioned in the Introduction, accu- 
rate computation of the GWs from EMRIs and the or- 
bital evolution is a complex and computationally inten- 
sive task. The most promising approach probably is the 
coupled integration of the compact object dynamics and 
GW emission taken in [S]. Alternatively, one could have 
a separate evolution of the orbital motion using self force 
computed across various geodesic orbits and employ os- 
culating elements approach [8j [14] . The waveform at in- 
finity could be obtained from the Teukolsky equations 
|15) in time or in frequency domain |161 117j . 

The above methods are computationally expensive and 
several approximations were suggested. Less accurate 
but still quite reliable are Numerical Kludge (NK) wave- 
forms: original NK [TT] and extended/improved NK 
called "Chimera" [18j EH] • Those methods combine accu- 
rate prescription for the orbital evolution with approxi- 
mate (Post-Newtonian) waveform generation formalism. 

The less precise model, which captures all relevant 
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physics of EMRIs (orbital precession, three orbital fre- 
quencies) was suggested in pQ , so called Analytic Kludge 
waveform. These waveforms are very fast to generate, 
and even though they cannot be used for searching for 
actual GW signals, they are used to develop data analysis 
algorithms and to evaluate their performance [U El [20] . 

In this work, we used NK waveform. In the original 
paper, the waveform was generated in the time do- 
main, we have reimplemented it in the frequency domain 
following suggestions of S. Drasco who did it first (private 
communications). However, we still take into account the 
time dependence of only three orbital constants under 
radiation reaction: energy, azimuthal component of the 
orbital angular momentum and Carter constant. Under 
the self force in osculating element approach we should 
evolve also other three constants (defining initial posi- 
tion of the compact object) due to conservative part of 
the self- force [5J [2] . This does not affect our search re- 
sults, since PW is model independent, however, we have 
to use the same model (as in the simulated data) for map- 
ping the phenomenological parameters onto the physical 
parameters of the binary. Mismatch in the models would 
result in the bias which we want to avoid. 

Finally we want to avoid using in this work the Ana- 
lytic Kludge model, because it predicts somewhat simpli- 
fied (detectable) harmonic content of the waveform. The 
NK waveforms for generic orbits were compared against 
waveforms based on the Teukolsky equation and they 
show quite good agreement. We believe that NK de- 
viates from the true EMRI signal in the phase but not so 
much in the number and strength of harmonics. There- 
fore we use NK model as a representation of the EMRI 
signal throughout this paper. 



III. EMRI PHENOMENOLOGICAL 
WAVEFORM FAMILY 

There are several algorithms which have been proven to 
be successful in detecting EMRIs in the simulated LISA 
data [13] [20] HI]. However, those algorithms partially 
utilize the features of AK waveform which was used in 
the simulation of the data and in the data analysis. As 
explained in Section [IT] we want to avoid it by building 
a generic phenomenological template family. 



A. Phenomenological waveform in the source frame 

The model we want to propose is based on the following 
assumptions about GW signals from EMRIs: 

1. The orbital motion can be effectively described by 
six slowly changing quantities. Explicitly, three 
time-dependent initial phases are governed by the 
conservative part of the self force; three fundamen- 
tal time-dependent frequencies are governed by the 
radiative part of the self force. 
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FIG. 1. The time- frequency plot of a typical EMRI signal 
without noise. There are 30 dominant harmonics in total. 



2. The waveform is represented by harmonics of three 
frequencies (phenomenologically, these frequencies 
are the summation of the fundamental orbital fre- 
quencies and the evolution of the initial phases) 
with slowly changing intrinsic amplitude: 



h(t) = ^ hlmn(t) 

l : m,n 

= Re j J2 A lmn {ty^A 

\l,m,n J 

= Re ( A lmn (t)e^+"^+^) j , (i) 

where $ r , $g, are the phase evolutions corresponding 
to the three fundamental motions. Here we omitted the 
tensorial spatial indices for simplicity. 

The first assumption basically expresses that the mo- 
tion is described by a slow drift from one geodesic to 
another. The initial phases correspond to the initial po- 
sition of a compact object on a given geodesic and the 
orbital frequencies are functions of the energy, azimuthal 
component of the orbital momentum and Carter con- 
stant. The slow drift ensures that phases <fri m n are slowly 
varying functions of time. 

Fig. [T] shows the time-frequency plot of a typical EMRI 
signal. There are 30 clearly separated frequency tracks in 
the noiseless plot, which display the dominant harmonics. 
It can also be seen that the frequencies of harmonics are 
smooth and vary slowly. It is generally true that both 
amplitude and the phase are slowly varying functions of 
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time, thus we can safely make the Taylor expansion: 

$ r (t) = $ r (t ) + w r (t )(t - to) + ^6j r (t - t ) 2 + ... 

= $ r (t ) + 27r/ r (*o)(* - *o) + w/r(* - to) 2 + . . . , 

(2) 

* fl (t) = $e(t ) + we(t )(t - to) + rWfl(* - to) 2 + • • ■ 

= $e(t ) + 2nfg(t Q )(t - t ) + nf e {t - t ) 2 + ..., 

(3) 

* v (t) = * y (to) + w^toX* - t ) + ^(t - to) 2 + ■ • ■ 

= $„(t ) + 27r/ v (to)(t - to) + ?r/ v (* - t ) 2 + ■ • ■ , 

(4) 

Almn(t) = A[ mn (to) + A lmn (t )(t - t ) + (5) 

Since the amplitudes Ai mn are even smoother than the 
phase over extended period of time, and because the de- 
tection techniques are more sensitive to mismatch in the 
phase than in the amplitude, we can neglect the time 
evolution in the amplitudes and treat all of them as con- 
stant. It is a very good assumption over three months 
of the simulated data which we analyze in this paper. 
As for the phase expansion, we calculate the so-called 
fitting factor (FF) for the different orders of polynomial 
approximations of the phase to check the fidelity of the 
PW. Numerical results show that the Taylor expansion 
for three months data, up to / order, gives the FF around 
0.9, and up to / order the FF is larger than 0.999. So it 
is sufficient to expand the phase to / order. This is the 
phenomenological waveform family which we propose to 
analyze an EMRI signal. To summarize, the phenomeno- 
logical waveform is a summation of individual harmonics 
with constant (or linear) amplitudes and polynomial (in 
time) phases. 

B. Prom the source frame to the LISA frame 

First we will express the GW wavefrom in the solar sys- 
tem barycenter frame and then translate it to the frame 
attached to LISA (or a LISA-like space based observa- 
tory). In the source frame, an arbitrary gravitational 
wave (GW) signal in the TT gauge can be written in the 
following form: 

h{t) = hl{t)e + + hi(t)e x (6) 

where the superscript 'S' denotes the source frame. Since 
the LISA constellation is orbiting the sun, it is convenient 
to express the GW signal in the solar system barycenter 
(SSB) frame. 

h(t) = h+(t)e + + h x (t)e x (7) 
e+ = 9 s ® 9 s - j> s <g) 4> s (8) 
e x = 9 s <g> <j> s + 4> s <g> 9 s (9) 



where (9 s ,4> s ) denotes the direction of the GW source 
in the SSB frame, 9 s , 4> s are the unit vectors along lon- 
gitudinal and latitudinal directions. The principal polar- 
ization vectors attached to the solar system barycenter 
frame, 6 , (j) S are connected to the principal polarization 
vectors in the source frame via rotation angle ip (since 
they lie in the same plane orthogonal to the GW propa- 
gation direction). The polarization components h + and 
h x are transformed under this rotation according to 

h+ = hi cos(2V>) + h% sin(2V>) (10) 
h x = -hi sin(2-0) + h x cos(2V>). (11) 

Now we will add LISA response. LISA measures the 
Doppler shift of the inter-spacecraft lasers induced by a 
gravitational wave. The single-link full response to this 
frequency shift can be derived with the help of three 
Killing vectors 22 1. However, this single-link signal is 
orders of magnitude smaller than the dominating laser 
frequency noise. Thus, we need to use the so-called Time- 
Delay-Interferometry (TDI) variables [23], which cancel 
the laser noise through the recombination of the artifi- 
cially delayed single-link signals. In the low frequency 
limit, the two orthogonal TDI (noise independent) vari- 
ables of Michelson type can be expressed as |5S] 

hj(t) = [SLi(t)-SL2(t)]/L 

= HO ■ Dj (12) 

hjj(t) = ^[5L 1 {t) + 8L 2 {t)-28L z {t)]/L 

- h(C) : D n (13) 

where L stands for the average arm length. The retarded 
time C(t) = t — k • x/c defines the wavefront, where k is 
the GW propagation direction. The two detector ten- 
sors are defined as Dj = \{n\ ® hi — n 2 ® ^2), Dn = 
■^j^(hi®hi+fi2®h2 — 2h 3 ®h 3 ), where 711,712,713 denote 
the unit vectors along each arm of LISA. Here we assume 
LISA-like setup which has six links (three arms). Even 
though the EMRI signal could reach quite high frequen- 
cies and require full response, we adopt the low-frequency 
approximation for our exercises. This does not restrict 
ability of our analysis as long as the simulated signal and 
the search template use the same response function. 

C. Data analysis with phenomenological waveform. 

We start with a brief overview of our notations and 
basics of data analysis. We denote the Fourier transform 
of a time series a(t) by a(f) and adopt the following 
convention 




We assume that the detector is characterized by a Gaus- 
sian, stationary noise n(t) and its two-sided noise power 
spectral density is defined as h*(f')h(f) = S n (f)S(f — 
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/'), where the over bar denotes the ensemble average. 
With this power spectral density, it is conventional to 
define an inner product of two time series a(t),b(t) as 
follows 



< alb >= 



, S n (f) 

The signal-to-noise ratio is defined as 



df. 



SNR 2 =< h\h >-- 



Sn(f) 



df. 



(15) 



(16) 



where h is the GW signal. Let us denote the proba- 
bility of a gravitational wave signal h(9) being present 
in the data s(t) by P(s\h(8)), where 9 is the set of pa- 
rameters that characterizes the gravitational wave signal. 
Similarly, the probability of no gravitational wave signal 
present in the data s is denoted by P(s\0). Likelihood 
ratio A(<?) is the ratio between these two probabilities 



Me) 



P(s\h(9)) 
P(s\0) 

_ e <s\h(9)>-±<h(9)\h(9)> 



(17) 



It is conventional to consider rather logarithm of the like- 
lihood ratio as a detection statistic: L(9) = log A(9) =< 
s\h(9) > -1 < h(9)\h(9) >. This is the quantity wc 
want to maximize over the parameter set 9. 

The likelihood ratio could be further simplified if we 
use PW. A single harmonic with polynomial phase up to 
/ order in the source frame takes the following form 



h(i) = A + cos($(t) + $ )e+ + A x sin($(t) + $ )e> 



$(t) = 2nf(t - t ) + nf(t - t Q ) 2 + 
£/(*-*d) 3 + ^7(*-fo) 4 , 



(18) 



(19) 



where we have omitted harmonic indices l,m,n. After 
simple algebra, LISA's response to this single harmonic 
GW signal without noise can be put in a simple form 

hj(t) = A%(t), hu(t) = A%\t) (20) 

where we follow summation convention over repeated in- 
dices, and (i = 1,2,3,4. The four amplitude parame- 
ters A^ depend only on (A + , A x , $ , ip), which are usu- 
ally called extrinsic parameters, while hAt),h* (t) are 

functions of (9 s , s , /, /, /, /), which are usually called 
intrinsic parameters. From now on, we denote the in- 
trinsic parameters by 9. The extrinsic parameters (be- 
ing constants in our approximation) can be maximized 
over analytically [551 HZ] > which we will show explicitly 
below. We denote the measured data with noise corre- 
sponding to hj(t),hu(t) by Si(t), Sn(t). Since the joint 
probability of a GW signal present in both sj and S/j is 
just the product of the individual probabilities, the joint 



log likelihood is just the summation of the individual log 
likelihoods 

L{9,A») = < 8i\hi(B) >-\< hi(9)\h!(9) > 

+ < s n \h u {0) >-\< h n (9)\h H {9) £21) 



Substituting ( 20 1 into this expression we arrive at 



L(9,A^=A^s I fl (0)-^M^(d)A v 

+ A»s»{9)- l -A»Mll{9)A\ (22) 

where we have used the following conventions: =< 

siK >> s " =< s^K 1 >> M U =< h lM >> M H =< 

hlf\hj, >■ We can maximize the log- likelihood over ex- 
trinsic parameters by solving 

dL(0, A^) 



dA^ 



(4 + *?) -(Mf u , + M£)A V = 0,(23) 



which is straightforward to find A^ = [{M + 
M ) _1 ] MI ' (sj, + si 1 ). The log-likelihood maximized over 
the extrinsic parameters is called F-statistic: 

F<9) = max L(0, A") 

= \(4 + 4')i( Ml + M")- l r(4 + 4 J ).(24) 

Its expectation value is connected to the SNR in the fol- 
lowing way 



E[F(9)} = ^SNR 2 + 2. 



(25) 



Since h(9) is narrow band signal, the inner product can 
be written in the following form 



<a\b> = J 



OO 

1 

Sn(fo) JQ 



a*(f)Hf) 

Sn(h) 
T 



df 

a(t)b(t)dt 



(26) 



where T is the observation time, /o is the middle fre- 
quency of h(6). The inner product is a function of T, 
and so is F-statistic. By varying T from to the to- 
tal observation time, we define a cumulative F-statistic 
F(T, 9). The cumulative F-statistic for 30 dominant har- 
monics without detector noise is plotted in Fig. [2] The 
case with the simulated detector noise is shown in Fig. |3j 
the total SNR of the signal in this case is SNR = 50. 
Those are two data sets which we will analyze in the 
next section. 

The cumulative F-statistic provides much more infor- 
mation than F-statistic. Actually, if 9* is the true pa- 
rameter set of the signal, one can argue that 



E 



dF(T, 9,) 
dT 



cx h 2 (T)e(T), 



(27) 
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Cumulative F-statistic of 30 harmonics Cumulative F-statistic of 30 harmonics, SNR = 50 




FIG. 2. The cumulative F-statistic of 30 dominant harmonics 
with true parameters without noise. Since there is no noise, 
the F-statistic is not normalized. 



FIG. 3. The cumulative F-statistic of 30 dominant harmon- 
ics with true parameters and detector noise. Note that the 
F-statistic is converted to SNR in the figure. The strong 
harmonics are cumulating gradually with local spikes. The 
low-SNR harmonics behave similar to noise, hence made un- 
detectable. 



where £(T) = y £+(T) + £| (T) is the geometrical mean 

of the antenna pattern functions for two polarizations. 

When there is no detector noise, dF ^jf^ = E dF ^ e ^ 

is nonnegative. Thus, E[F(T,0*)] is always increasing 
over the entire time span when the GW signal is present, 
as can be seen in Fig. [2] It is not necessarily so in presence 
of the noise and during analysis of the data. There are 
three types of oscillations on the cumulative F-statistic 
curve F(T,0). (i). The (non-negative) oscillation due to 
the oscillatory nature of the gravitational wave signal. It 
is at twice the GW frequency, which makes it hard to see 
m Fig. [2} (ii). In reality, we do not know the exact true 
parameters of the GW signal. That means, in most cases, 
the parameter set 6 we try differs from the true param- 
eter set 0*. This introduce beat-notes to F(T,0). This 
kind of oscillation happens at beat-note frequency, which 
is much lower than the GW frequency itself, (iii). The 
third type of oscillation is due the noise. The presence 
of the noise makes the cumulative F-statistic uneven, see 
Fig. [3j Comparing to the former two types, this kind of 
oscillation is irregular; it oscillates at all frequencies and 
could cause temporary (for a short time) decrease in the 
cumulative F-statistic. 



We have found that over three months of simulated 
data we can consider all harmonics as being completely 
independent with virtually zero overlap between them, 

< hmn\h'm'n' >= $W 5 mm > $nn> ■ The total F-statistic is 

therefore a sum of F-statistics from each harmonic. In 
the next section we describe the search where we use Eq. 
(24) as a detection statistic, and we will use cumulative 



F-statistic later on to analyze our findings. 



IV. SEARCH WITH THE 
PHENOMENOLOGICAL WAVEFORM 

In this section, we use the PW as described above 
together with the introduced detection statistic. We 
will use two 3 month worth simulated data sets: with 
and without noise. We use the same GW signal (based 
on NK model) in both cases. The total SNR of the 
source in the noisy case is 50. We have taken the fol- 
lowing parameters for the EMRI : the mass of the MBH 
M = 10 6 M Q , the mass of the compact object (stellar 
mass BH) m = 1QM@, the initial orbital eccentricity 
e = 0.4, the semi-latus rectum p = 8M, the inclina- 
tion angle i = ir/9, the spin of the MBH a = 0.9M, the 
sky position of the source (8 s ,(f) S ) = (n/i, 7r/4), the po- 
larization angle tp = 0. In our analysis we assume that 
the sky location is known. Our primary goal here is to 
recover the intrinsic parameters of the source. 

The noiseless case is used to avoid any possible bias 
in the final result due to stochastic nature of the noise, 
and assess possible restrictions of our search technique 
and PW family. Next, we apply the same search method 
to the same GW signal buried in the noise, which would 
justify its effectiveness in practice. 

Here, we describe the search for individual harmonics 
with Markov chain Monte Carlo (MCMC) method. For 
completeness and future references we give a brief intro- 
duction to MCMC. Like a standard Monte Carlo integra- 
tion, MCMC is a random sampling method. It is noth- 
ing but Monte Carlo integration with a Markov chain. 
By properly constructing a Markov chain, MCMC can 
draw samples from the searching parameter space more 
efficiently. Among all the methodologies of constructing 
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a Markov chain, the Metropolis-Hastings scheme would 
be the most general one. The main idea of Metropolis- 
Hastings algorithm is to cleverly construct a Markov 
chain that satisfy the detailed balance equation, so that 
the sampling distribution will converge to the likelihood 
surface we want to estimate. If the shape of the likelihood 
surface is known, the parameter set that corresponds to 
the maximum likelihood is automatically known. Thus, 
MCMC is also widely used as a stochastic optimization 
tool in GW data analysis (we refer the reader to a very 
nice overview and discussion on Bayesian methods in |28j , 
see also references therein). 

If the likelihood surface is multimodal(i.e. contains 
large number of separated local maxima) then simple 
version of the MCMC finds a maximum and does not 
move off it to explore larger parameter space. Many ways 
around this problem were suggested but we will not use 
any of them here (besides simulated annealing which we 
will discuss a later). As we will see, a simple Metropolis- 
Hastings algorithm is sufficient. The likelihood surface 
of an EMRI signal is very rich in "wall" and "needle" 
like structures, which make it very hard to find a global 
maximum. We are interested in detecting as many local 
maxima as possible. Therefore we run multiple indepen- 
dent chains and harvest the results after they converge to 
various maxima of the likelihood surface. If we are lucky, 
the global maximum could be among multiple maxima 
we have found. 

To understand the Metropolis-Hastings algorithm, 
first consider a stochastic process denoted by {9k\k = 
0, 1, 2...} which belongs to the parameter space B in W 1 . 
Here we defined k as a set of parameters at step fc, which 
can also be viewed as a point in the parameter space B. 
If there exists a transition probability P(0 k +i\0k) de- 
pending only on the current point k for the stochastic 
process to be in state 8k+i, we call this stochastic pro- 
cess {0 k \k = 0, 1, 2...} a Markov chain with a transition 
probability P{8 k+ i\0 k ). In a Bayesian viewpoint, we can 
take this transition probability as conditional probability 
and immediately see that 



/ p(e k+1 \d k )d0 k+1 = i . 

■IB 



(28) 



A Markov chain satisfying the detailed balance equation 



A(0 k )P(0 k+1 \0 k ) = A(0 fc+1 )P(0 fe |0. 



k+l) 



(29) 



will (up to some relatively weak conditions) be equivalent 
to the samples from the distribution A(0) after a certain 
initial period (often called burn-in stage) . We can easily 
estimate the distribution A(0) with the Markov chain 
samples and hence the most probable parameter set 6 
for given observed data s, where 



A(0|s) 



max A(0|s) 
e 



(30) 



is usually called the maximum likelihood estimator. 
By virtue of Metropolis-Hastings algorithm, we can con- 
struct a Markov chain that satisfies the detailed balance 
equation and make use of the corresponding property to 
estimate our template parameters 0. To do this, we ran- 
domly choose a parameter set 0q in the parameter space 
as the starting point. Then one can pick a proposal dis- 
tribution q(9 k+ i\6 k ) (as long as there is no forbidden re- 
gion in the prescribed parameter space to the point k +i) 
and sample a candidate point 9 k+ i from this distribution. 
Then we calculate the acceptance probability defined by 
the following formula 



a(8 kl k+1 ) = min 1, 



A(0 k+1 )q(0 k \0, 



fc+l. 



A(0 k )q(0 k+1 \0 k ) 



(31) 



By accepting the point k +i according to the above prob- 
ability, we have, in fact, succeeded to construct a transi- 
tion probability, 



P(0 k +l\0 k ) = q(0 k+ i\0 k )a(9 k ,0 k+1 ) 



(32) 



It is easy to see that the Markov chain generated by the 
above transition probability satisfies the detailed balance 
equation: 



A(0 k )P(0 k+1 \0 k ) = mm{A(0 k )q(0 k+1 \0 k ),A{0 k+1 )q{0 k \0 k+1 )) 
= min (A(0 k+1 )q(0 k \0 k+1 ),A{0 k )q{0 k+1 \0 k )) 
= A(0 k+1 )P(6 k \0 k+1 ). (33) 

Thus, such a Markov chain will eventually serve as a 
succession of samples from A(0). The best performance 
is achieved if the proposal probability q(0 k +i\0 k ) resem- 
bles the target distribution A(0) over the entire param- 
eter space. Without prior knowledge about the kind of 
probability distribution around the true parameter loca- 
tion, it is natural to choose it as a multivariate normal 
distribution centered at the present point k with covari- 
ance matrix C, 



q(0 k +l\0k) 



v /(27r) Jv det[C] 



exp 



,T C 1 (0 k +i — 0k) 
(34) 



where N denotes the dimension of the parameter space 
and det[C] the determinant of the covariance matrix C. 
The likelihood surface has usually multimodal (multiple 
local maxima) structure, and, therefore, a single multi- 
variate normal distribution cannot describe the proba- 
bility density over the entire template space but only a 
very small region around the local maximum. Since the 
probability distribution at the local maximum is usually 
very sharp, a Markov chain easily gets trapped there for 
many steps. To avoid insignificant maxima we use the 
so-called annealing scheme, originating from simulated 
annealing. We adopt two types of annealing techniques: 



(i). we introduce a temperature 71 to the acceptance 



rate a [equation (31)] so as to have a larger possibility 
to accept the proposal point in the beginning. By com- 
bining equations (IT), (24), (31), (34), the acceptance 



probability is now written as 



a(8 k ,9 kA 



min ( 1 , e 



,[F(e k+ i)-F(0 k )]/Ti 



(35) 



where the temperature 71 = 71 (fc) is a function of the 
step index k, it starts from some relatively large num- 
ber and gradually decays to unity, (ii). We introduce 
a second temperature 7i to the proposal distribution 
q(0k+i\9k)- The covariance matrix C is replaced by 
C x T%- Same as 71, 72 is also a function of the step 
index k, decaying gradually to unity. Hence, the chain 
take large steps in the beginning and explores large vol- 
ume in the parameter space. Explicitly, we choose 71 and 
72 both as a linear function of k with negative slope. 
Let us summarize the algorithm: 

1. k = 0. Choose a random parameter set 8q as the 
starting point and calculate the F-statistic F(0q). 

2. k — > k + 1. Calculate the temperature 71(fc),72(fc)- 

3. Generate the next candidate parameter set C from 
the proposal distribution with modified covariance 
C xT 2 . 

4. Calculate the F-statistic of the new parameter set 
F(8 C ). 

5. Calculate the acceptance probability a(0k,0 c ) = 
min(l, e [ F ^)-F(e fc )]/7i). 

6. Draw a random number u from unity distribution 
U(0, 1). If u < a, accept the candidate parame- 
ter set 6k+\ = Cl else, stay at the current point 
9k+i = Sk- 
in the search we have used a diagonal form of 
the covariance matrix in the gaussian proposal dis- 
tribution (34), with the following elements: C = 
[diag(l(T 4 , 1CT 12 , 1(T 20 , 10~ 28 )] 2 corresponding to the 
parameter set {/, /, /, /}. And 72 used to scale the co- 
variance matrix decays linearly with the number of mem- 
bers in the chain from 1 to 5 x 10~ 4 . We have found 
that the use of the actual Fisher information matrix as 
C did not improve significantly the search results. We 
run about 50 chains on both noiseless data and noisy 
data. All the parameter sets that generate an SNR larger 
than a certain threshold (we have used SNR > 4.5) are 
recorded. Notice that there are possibly many such qual- 
ified parameter sets in a single chain. Thus, we have 
hundreds to thousands of qualified parameter sets or lo- 
cal maxima. These local maxima contain information 
about the signal. We will analyze these local maxima in 
the next section. 



V. ANALYSIS OF THE SEARCH RESULTS AND 
MAPPING TO THE PHYSICAL PARAMETERS 

In this section we will explain how we use the results 
of MCMC search described in the previous section and 
reconstruct harmonics of the GW signal. Furthermore, 
we use the model of EMRI (NK) to estimate the physical 
parameters of the system. 



A. Clustering algorithms 

In this subsection we extract information from the lo- 
cal maxima detected by MCMC search. We first focus on 
the noiseless data to explain the algorithm, then modify 
it a bit and apply it to the noisy data. Since this work is 
the first of a series of papers, the main task here is to es- 
tablish the framework and justify the method. Hence, as 
mentioned above, we have assumed that the sky position 
of the source is known and concentrate on the intrin- 
sic parameters only. This will save us some time, yet 
maintain all the main features of the problem. As a re- 
sult, each local maximum is characterized only by the 
frequency and its derivatives (/,/,/, /)• 

Let us look at one example to understand how we 
extract the information about the source from the de- 
tected local maxima. We take a particular solution of 
MCMC search and for each harmonic of PW we can 
compute cumulative F-statistic according to the prescrip- 
tion given in Section [ill C| We concentrate only on those 
harmonics which give significant contribution to the to- 
tal F-statistic. If the harmonics of PW match perfectly 
the harmonics of a signal we should observe something 
similar to Fig. [3j however it is rare when we detect a 
full harmonic (only sometimes for the strongest). More 
frequently, we detect a part of a harmonic (frequency 
and derivatives close to true but not exact) or even sev- 
eral harmonics at different instances of time as shown in 
Fig. [4] The black and green curves are two strong har- 
monics of a signal (black being stronger), and the blue 
is a harmonic of PW. In the pink regions, our template 
matches for a short period of time the frequency of a 
signal (two distinct harmonics at two instances). The 
corresponding cumulative F-statistic is shown in Fig. [5] 
There are two positive jumps in the accumulation of the 
F-statistic which correspond to two instances of inter- 
section. Therefore, we can conclude that the positive 
slope in the cumulative F-statistic (if it happens over a 
significant duration) corresponds to the part of the fre- 
quency and time where a harmonic of PW matches (at 
least partially) some harmonics of a signal. We collect 
such events of matching and display them on the time- 
frequency plane, resembling the mosaic of a true signal. 

The violent oscillation in Fig. [5] is one of the three 
types of oscillations on the cumulative F-statistic curve 
mentioned in the previous section. In fact, it is the beat 
note between the true harmonics and the local maximum. 
Observe that the beat notes happen at relatively higher 
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FIG. 4. Time-frequency plot of harmonics. The black and 
green tracks are two strong harmonics of the EMRI signal 
(black being stronger). The blue track corresponds to a har- 
monic of PW that accumulates a significant F-statistic. It 
intersects the true harmonics at the pink segments, those cor- 
respond to times of increase of F-statistic, see Fig. [5] [6] 



frequency, while the increasing slopes (where the local 
maximum matches the frequencies of the true harmon- 
ics) have relative low frequency. Thus, we design a third- 
order Butterworth low pass filter to get rid of the beat 
notes. After the low-pass filter, the cumulative F-statistic 
has only few extrema, as shown in Fig. [6j After clearing 
up the cumulative F-statistic, we apply two criteria for 
identifying a significant F-statistic accumulation: (i) the 
slope must be larger than certain threshold; (ii) the accu- 
mulation time must be over longer than certain period. 
As it is seen by eye tuning those two parameters should be 
sufficient to get the right parts of cumulative F-statistic. 
In our search we have made the following choice for those 
parameters. In the case of noiseless data, we require the 
slope to be larger than one-tenth of the largest slope of 
the cumulative F-statistic of that trial harmonic, and the 
cumulative time (over which we observe steep positive 
slope) to be longer than three days. 

We plot all recovered patches on the time-frequency 
plane in Fig. [7| where we can identify by eye 13 strong 
harmonics. Although the weaker harmonics are lost, the 
strong ones retain enough information about the EMRI 
system evolution, hence allowing us to recover the physi- 
cal parameters we are interested in. Zooming at a specific 
harmonic in time and frequency, one will see that there 
are many patches from different results and at each in- 
stant we observe a finite spread in the frequencies for 
a given harmonic. This is due to various solutions from 
MCMC search matched a given harmonic of a signal with 
different precision. However, we expect that the distribu- 
tion of found frequencies at each instant of time will be 
centered at the true frequency of the signal's harmonic. 
As an example, we show distribution of found frequen- 



FIG. 5. Unfiltered cumulative F-statistic corresponding to 
the PW harmonic and data given in Fig. [4] The F-statistic 
labeled on the vertical axis has only relative meaning, since 
we work with the noiseless data. The green and red squares 
mark the extremes of the curve, thus distinguishing between 
the increasing and the decreasing slopes. The large number of 
the extremes is due to the beating between the true harmonics 
and the trial harmonic. 
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FIG. 6. Filtered cumulative F-statistic corresponding to 
the situation depicted in Fig. [4] It is similar to Fig. [5] but 
after applying the low pass filter to remove the beatings (high 
frequency oscillations). 



cies at a particular instance of time for two harmonics 
in Fig. [8j In that plot we show the histogram of de- 
tected frequencies at that time in blue and Gaussian fit 
as smooth green curves. This is to be compared with fre- 
quencies of two harmonics of a signal at the same time in 
red. As mentioned above, different solutions of MCMC 
search vary in precision of matching the signal at different 
instances, and we can use accumulation time as a mea- 
sure of goodness of match of a signal by a given solution. 
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FIG. 7. Time-frequency plot of all patches corresponding 
to strong accumulation of F-statistic. We can identify parts 
of frequency tracks of 13 EMRI harmonics. Each track in 
this plot has a finite width coming from different solutions of 
MCMC search which have different precision of matching the 
signal. 



FIG. 8. Zoom at two harmonics at a specific instance of time. 
The red stems denote the frequencies of the true harmonics of 
a signal, while the blue histogram shows the detected frequen- 
cies at this instant. The green curves display the Gaussian fit 
to the frequency data with re-scaled amplitudes. The verti- 
cal axis of pink points indicates the relative time over which 
we have observed strong accumulation of F-statistic for each 
solution. 



The relative accumulation time of different solutions are 
shown as pink points in Fig. [H] First, one can see that 
Gaussian fit lies on the top of the true frequency, and 
second, that the distribution of pink points is similar to 
the blue histogram, so either can be taken to character- 
ize the found harmonics of a signal. Similarly, we can 
do at each instance of time for all found tracks in the 
time-frequency plane. For the noiseless search we picked 
uniformly 10 instances and made a Gaussian fit around 
each harmonic. We identify the mean of the Gaussian fit 
as the most likely frequency of a signal's harmonics at 
that instance and we identify the spread (standard devi- 
ation) of a distribution as an error in our evaluation of a 
frequency. The result of this clustering is given in Fig. [9] 

In the case of data with the detector noise, the ba- 
sics and the strategy are roughly the same as in the 
noiseless case with minor modifications. In the begin- 
ning, we record the local maxima with SNR greater than 
4.5. Next, we select the significant increasing slopes of 
the cumulative F-statistic with three requirements: (i) 
the maximum F-statistic along the cumulative F-statistic 
curve is larger than 50, (ii) the minimum slope of the 
significant increasing segment is larger than 4 x 10 s , 
(iii). the monotonic increasing duration is longer than 
about a week. Those conditions are more stringent than 
for the noiseless case and eliminated several found weak 
harmonics of EMRI signal. However, at the same time 
they significantly reduce the false events (and that is 
what we want). From this selection, we identify 5 strong 
harmonics in the noisy case. After that the procedure is 
similar to the noiseless case. 
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FIG. 9. Gaussian fit to the detected frequencies at ten 
instants. The red points represent the mean of a Gaussian fit 
as shown in Fig. [§] for each harmonic at ten instants . The 
blue error bars show the la uncertainties of the Gaussian fits. 
Note the tiny error bars are along the frequency dimension 
which indicates that the MCMC search localizes quite well 
frequencies of the EMRI's harmonics. 



B. Search for physical parameters 

Now we are in a position to recover the physical pa- 
rameters of the binary system. First, we need to adopt 
the model for the orbital evolution, and here we used the 
same model as used in the simulation of the data sets. 
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In the noiseless case the only reason for the deviation of 
recovered parameters from the true values is due to inac- 
curate identification of the tracks in the time-frequency 
plane or due to ambiguity in solving the inverse problem 
(mapping harmonic tracks onto the physical parameters, 
m/M, a, e, i,p/M). We have performed the search on the 
time- frequency plane similar in spirit to 29J. We have 
used weighted chi-square test 




between signal tracks (for different parameters) and 
recovered tracks (Fig. [9| . We have used particle swarm 
optimization (PSO) and genetic algorithm (GA) as two 
independent search methods to test the robustness of our 
result. We start with describing the PSO method, and 
then give brief overview of GA. 

Particle swarm optimization (PSO) is a stochastic opti- 
mization method introduced by Kennedy and Eberhardt 
in 1995 [30]. In gravitational wave data analysis, PSO 
was first applied to a binary inspiral signal [31 j . In this 
section, we briefly describe the algorithm, while further 
details can be found in the references [3jTl |3T] • 
The goal of PSO is to find the global minimum/maximum 
(here we minimize the chi-square test) of a parameterized 
functional k(6) and the corresponding parameter set 0*, 
where 6 stands for an arbitrary parameter set in W 1 . The 
idea is to evaluate simultaneously at different pa- 

rameter sets 6i, i — 1, 2, treating them as particles in 
the parameter space, and evolve them according to cer- 
tain dynamics until the stable solution is reached. Let us 
denote the z-th particle out of a swarm of N p particles 
during fc-th iteration in the search by 9i[k]. Its position 
in the parameter space in the next iteration is determined 
by its velocity in the current iteration vi [k] , 

9 i [k + l]=9 i [k}+v i [k}. (36) 

Usually, the particles start with randomly chosen posi- 
tions Oi[l) and velocities Vi[l]. Up to fc-th iteration, we 
denote the i-th particle's best location by Q\ [k], in the 
sense that 

K (0f[fc])=min«(0 4 [j]). (37) 

The global best location Of [k] up to the fc-th iteration is 
defined by 

K (6 9 [k}) = min/t(0f[fc]). (38) 

i 

Note that particle best locations and the global best lo- 
cation are the best parameters respectively found by the 
individual particles and the whole swarm in the entirely 
history of the search up to the fc-th iteration. They 
are updated only when a better parameter set is found. 
These best locations contain a lot of information about 



the functional k(9), so they are used to guide the parti- 
cle's motion in the future. Explicitly, the velocities are 
updated with the following equation 

Vi[k + 1] - wvi[k] + cixi(0f [k] - e t [k]) + 

c 2 X 2 (0 9 [fc]-0,:[fc]), (39) 

where w is called the inertia weight, C\, C2 are called 
the acceleration constants (we take them to be the same 
as in |3T]) and xij Xi are random numbers drawn from 
U(0, 1). We run PSO search several times until the 
return result is confirmed by several searches. 

The second search method is called Genetic Algorithm 
(GA) and there we evolve a number of parameter sets 
(points in the parameter space R n ). Each parameter set 
6i is called an organism, individual parameters are called 
the genes of this organism and the set of organism at fc-th 
search iteration step is called fc-th generation. We evolve 
generations according to the prescribed rules called "par- 
ents selection" , "breading" and "mutation" . The main 
idea of this optimization technique is to evolve colony of 
organisms toward the better fitness (which could be like- 
lihood ratio or, in our case, chi-square value) like in Dar- 
win's theory of natural selection. The strong organisms 
(with better fitness) participate more often in breading 
and therefore drag the colony toward the better values 
(lower) of chi-square. Mutation brings element of ran- 
domness in the search and occasional "positive" muta- 
tions help to avoid trapping around local minimum. For 
use of GA in GW data analysis we refer to [3H [33J and 
references therein. 

Let us give few more details specific to the implemen- 
tation used in here. We use x 2 value as a measure of 
fitness for each organism (smaller value is better). In 
each generation we use the roulette method with the se- 
lection probability proportional to the fitness of each or- 
ganism. For breeding we have used the one random point 
crossover rule. The probability mutation rate is mono- 
tonically decreasing function of the generation number: 
we have started with high probability of mutation to ex- 
plore a large part of the parameter space and decrease 
it gradually as organisms converge to a particular part 
of the parameter space. We have used "children" and 
"parents" sorted in the fitness to make a new generation: 
we use 50% of the best organisms. We automatically 
achieve the "elitism" in a way that the best % 2 value is 
never increasing from one generation to the next. 

We use the multi-step method to accelerate the search. 
In each step we evolve the colony for 500 generations 
as described above, but each new step uses the last 
generation of the previous step as the initial state. We 
have started evolution in the first step with completely 
random distribution of the organisms. The evolution 
of the colony at each step finishes with a very small 
mutation probability and with organisms confined to a 
quite small volume of the parameter space. The conse- 
quent search steps ensure that the found solution is a 
robust solution with respect to increase of the mutation 
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probability which disperses organisms forcing them to 
explore the parameter space for presence of a solution 
with better fitness. This helps to avoid being trapped in 
the local minima. The termination condition is the sta- 
bility of the best solution over several steps of the search. 

We have applied both those methods to fit the found 
tracks on the time frequency plane with the harmonics of 
EMRI signal. The search is done in 5 dimensional param- 
eter space with quite broad priors on (e,p/M, l, a, p = 
m/M), those are the eccentricity, the semi-latus rectum, 
the orbital inclination angle at the moment of beginning 
of observation, the spin of the MBH, and, the mass ratio 
between the stellar BH and the MBH. The total mass is 
not present here, we have kept it fixed to M — 10 6 M Q . 
For a given set of parameters, our search algorithm com- 
putes three fundamental orbital frequencies as functions 
of time, then a weighted chi-square goodness of fit test 
is preformed on harmonics of the signal. We use the 
means and standard deviations from the Gaussian fit as 
found point and its error in the time-frequency plane. 
The best fit corresponds to the lowest value of x 2 - We 
have used harmonics of the signal, which are expected 
to be strong over the large part of the parameter space, 
and have found this "harmonic table" by intensive montc 
carlo with NK models generated in the frequency domain. 
The index table has been truncated by choosing harmon- 
ics contributing (in total) 90% of the overlap with a total 
signal [j] 

The recovered parameters are given in the table [T] 



VI. SUMMARY 

In this paper we have introduced the phenomenological 
family of waveforms (PW) for detecting EMRI signals in 
the data from the LISA-like observatory. The template 
is constructed out of independent (over the time inter- 
val we have applied our analysis) harmonics of slowly 
evolving three orbital frequencies. We have neglected the 
amplitude evolution and presented the phase as a Tay- 
lor series up to the third derivative of frequency. Our 
analysis was restricted to the case of monotonically in- 
creasing frequencies. This condition will break only close 
to the plunge. The number of harmonics and range of 
indices were taken from the analysis of dominant har- 
monics of our model signal, though we have found at the 
end that the search only weakly depends on the number 
of used harmonics (only through the accumulated total 
SNR, which should be sufficient to claim detection). 



The total signal here to be a NK waveform with a large number 
of harmonics. We still truncate the number of harmonics used 
to build the signal: we stop if the inclusion of the next harmonic 
does not change overlap with the already built signal by more 
than 0.1%. 



Constructed phenomenological templates allows us to 
search for EMRI signals in a model independent way. 
This way we avoid complexity of accurate modeling the 
orbital evolution and gravitational waveform during the 
search. In addition PW cover also all possible small de- 
viations of the background spacetime from the Kerr solu- 
tion which would influence the signal's phase and could 
lead even to loss of the signal if the template assumes 
pure Kerr background geometry. 

We have used MCMC based search to find a large 
number of local maxima of the likelihood surface. We 
were not that lucky to find the global maximum. We 
have analyzed the found solutions by means of cumula- 
tive F-statistic over the time and identified the patches 
of the signal which were match by templates. As a re- 
sult, we have constructed a time- frequency map of (parts 
of) the signal's harmonics. Each track could be char- 
acterized by the best guess and the error bar at each 
instance of time (by fitting Gaussian profile to found fre- 
quencies around at that time each track). The next step 
is to assume a model for the binary orbital evolution, and 
check if the found time-frequency picture corresponds to 
the strongest harmonics of a signal. In other words, we 
want to find the physical parameters of the binary system 
which strong GW harmonics could leave the found im- 
print. We do that by conducting a search using particle 
swarm optimization techniques and, independently, ge- 
netic algorithm. We have used weighted chi-square good- 
ness of fit test to choose the best matching harmonics of 
the signal. We have assumed the same model as was used 
in the simulated data, and the recovered parameters are 
within 2% of the true values. 

We want to make few final remarks, (i) The found 
time-frequency tracks of the GW signal from EMRI did 
not assume any particular model. The mapping of these 
tracks to the physical parameters could be done in post 
processing using several models. We have chosen on pur- 
pose rather short (3 month) duration of the data. The 
search procedure could be repeated for each three months 
and then one can check consistency of a given model 
or further improve accuracy in the recovered parameters 
(if our model gives consistent parameters of the system 
across different data segments). This could be a power- 
ful method to search deviations from "Kerness" . (ii) In 
the mapping of the time-frequency tracks to the physical 
parameters of the binary, we have only weakly used infor- 
mation about the strength of each track/harmonic. We 
have found that the information stored in the frequency 
evolution is sufficient to recover parameters of EMRI. 
However, additional information about the strength of 
the recovered harmonics and harmonics of the modeled 
GW signal could give us additional confidence in the re- 
sult and/or distinguish between several solution, if am- 
biguity happens, (iii) Mapping from the found time- 
frequency tracks onto the physical parameters might turn 
out to be the most computationally intensive task. How- 
ever, one might use the information about the strength 
and a number of found harmonics to restrict a volume of 
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TABLE I. Recovered parameters of EMRi against actual parameters used in simulated data sets. 



description 


c(to) 


p(*o) 


fc(*o) a 


A* 




True parameters 


0.4 


8.0 


0.349 0.9 


10" 


5 


Recovered parameters (with noise) 


0.395 


8.029 


0.342 0.891 


9.79 x 


10" 6 


Recovered parameters (no noise) 


0.402 


7.991 


0.360 0.901 


1.002 x 


10" 5 



the searched parameter space. In addition, to perform 
mapping we require mainly the computation of the or- 
bital evolution, not the full waveform. However, it is then 
important to know which harmonics are the strongest for 
a given parameter set. (iv) In the future work we intend 
to include the sky location and the MBH mass into the 
search and investigate the possibility to differentiate be- 
tween different models of EMRIs based on the results of 
MCMC search with PW (as discussed in (i)). 
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